// Basu & Bundick (2018)
// "Uncertainty Shocks in a Model of Effective Demand: Reply"
// Alternative Preferences & Alternative Calibration
// Preference Specification is Equation (3) of the Reply with IES = 0.5

var a, c, de, deltau, deltauprime, expre, expre2, expvf1sigma, inv, k, mu, n,
	pe, pie, profit, q, r, rr, rrk, sdf, u, varexpre, vf, vola, w, y, z; // sdfplus; // // expvfconbundle1, expvfconbundleaux, expvfconbundle2, d2vfdc2, cn, dvfdc; 
    // dc,
    // cre, cre2, varcre, cre3, skewcre, cre4, kurtcre, marginalutility, cn, 
    // mpc1, marginalutility2;// , mpc2, mpc3, mpc4, mpc5, mpc6, mpc7, mpc8;
    // expQnum2_aux, expQnum_aux, expQnum2, expQnum, exsdf, risk_neutralR2, risk_neutralR, risk_neutral_variance, vrp;

varexo ea, evola, ez;

parameters 	alpha, beta, delta0, delta1, delta2, eta, phip, phik, 
			sigma, ies, frischlse, thetavf, thetamu, leverageratio, fixedcost,
			rhor, rhopie, rhoy,	rhoa, rhovola, volvola, rhoz, productionconstant, utilityconstant, 
            ass, css, dess, deltauss, deltauprimess, express, expre2ss, expvf1sigmass, invss, kss, nss, muss, 
            pess, piess, profitss, qss, rss, rrss, rrkss, sdfss, uss, varexpress, vfss, volass, volzss, wss, yss, zss;

phikestimate			= 4.1527;
rhovolaestimate 		= 0.8093;
volvolaestimate         = 0.003968; // 
volassestimate          = 0.0043648;
volzssestimate          = 0.018144;
rhoaestimate            = 0.9755;
rhozestimate            = 0.48083;

// phikestimate			= 3.9208;
// rhovolaestimate 		= 0.76889;
// volvolaestimate		= 0.0043735;
// volassestimate         = 0.0049905;
// volzssestimate         = 0.018792;
// rhoaestimate           = 0.98096;
// rhozestimate           = 0.48083;
   						
alpha           = 0.333;
beta            = 0.994;
delta0          = 0.025;
delta1          = (1/beta - 1 + delta0);
delta2          = 0.00031;
phip            = 240;
phik            = phikestimate;
sigma           = 100; // 100;   //40
ies             = 0.5; // 0.25
frischlse       = 2.0; // 2.5
thetavf         = (1 - sigma)/(1 - 1/ies);
thetamu         = 6; // 6.0;
leverageratio   = 0.90;

piess           = 1.005;
rss             = piess/beta;
rhor            = 0.000;
rhopie          = 1.5;
rhoy            = 0.2;

ass             = 1;
rhoa            = rhoaestimate;
volass          = volassestimate;
rhovola         = rhovolaestimate;
volvola         = volvolaestimate;

zss             = 1;
rhoz            = rhozestimate;
volzss          = volzssestimate;

// Steady State Calibration

yss         = 1.0;
rrkss       = 1/beta - 1 + delta0;
muss        = thetamu/(thetamu - 1);
uss			= 1;
fixedcost   = (muss - 1)*yss;
kss         = alpha*(yss + fixedcost)/(rrkss*muss);
invss       = delta0*kss;
css         = yss - invss;

nsssolution = fsolve(@(nssguess) 1+(-1+1/ies)/(1+((-1+alpha)*(-1+nssguess)*(fixedcost+yss))/(css*nssguess*muss)) -(frischlse/ies)*nssguess*(1-nssguess)^(-1),0.5,optimoptions('fsolve','Display','off'));
nss         = nsssolution;

wss                 = (1 - alpha)*(yss + fixedcost)/(nss*muss);
eta                 = (wss*(1 - nss)/css + 1)^(-1);
qss                 = 1;
productionconstant 	= (yss + fixedcost)/(kss^(alpha)*(nss)^(1 - alpha));
vfss                = 1;
utilityconstant     = vfss^((1 - sigma)/thetavf)*(1 - beta)/(css^(eta)*(1 - nss)^(1 - eta))^((1 - sigma)/thetavf);
dess                = css - wss*nss + leverageratio*kss*(beta - 1);
pess                = beta*(1-beta)^(-1)*dess;
profitss            = (muss - 1)*yss - fixedcost;
express             = 1/beta;
expre2ss            = 1/beta^(2);
varexpress          = 0;
expvf1sigmass       = vfss^(1-sigma);
qss                 = 1;
rrss                = 1/beta;
sdfss				= beta;
deltauss			= delta0;
deltauprimess		= delta1;

model;

y + fixedcost 			= productionconstant*(z*n)^(1 - alpha)*(u*k(-1))^(alpha);
c + leverageratio*k/rr 	= w*n  + de + leverageratio*k(-1);
w 					= ((1 - eta)/eta)*c/(1 - n);
vf 					= (utilityconstant*(a*c^(eta)*(1 - n)^(1 - eta))^((1 - sigma)/thetavf)  + beta*expvf1sigma^(1/thetavf))^(thetavf/(1 - sigma));
expvf1sigma 		= vf(+1)^(1 - sigma);
w*n 				= (1 - alpha)*(y + fixedcost)/mu;
rrk*u*k(-1) 		= alpha*(y + fixedcost)/mu;
q*deltauprime*u*k(-1) = alpha*(y + fixedcost)/mu;
// u = 1;
// k = kss;
k 			= ((1 - deltau) - (phik/2)*(inv/k(-1) - delta0)^(2))*k(-1) + inv;
deltau		= delta0 + delta1*(u-1) + (delta2/2)*(u-1)^(2);
deltauprime	= delta1 + delta2*(u-1);
 sdf = beta*((a/a(-1))*(c^(eta)*(1 - n)^(1 - eta))/(c(-1)^(eta)*(1 - n(-1))^(1 - eta)))^((1 - sigma)/thetavf)
			*(c(-1)/c)*(vf^(1 - sigma)/expvf1sigma(-1))^(1 - 1/thetavf);
//sdf = beta*(a/a(-1));
// sdf = beta*((a/a(-1))*(c^(eta)*(1 - n)^(1 - eta))/(c(-1)^(eta)*(1 - n(-1))^(1 - eta)))*(c(-1)/c);
// sdf = beta;//*((a/a(-1))*(c^(eta)*(1 - n)^(1 - eta))/(c^(eta)*(1 - n)^(1 - eta)))^((1 - 1)/thetavf)
		   //	*(c/c)*(vf^(1 - 1)/expvf1sigma(-1))^(1 - 1/1);

// sdfplus = rr*sdf(+1);
1 = rr*sdf(+1);
1 = r*sdf(+1)*(pie(+1))^(-1);
1 = sdf(+1)*(de(+1) + pe(+1))/pe;
log(r) = (1 - rhor)*(log(rss) + rhopie*log(pie/piess) + rhoy*log(y/y(-1))) + rhor*log(r(-1)); // 0.000000000001*rhoy*
de = y - w*n - inv - (phip/2)*(pie/piess - 1)^(2)*y - leverageratio*(k(-1) - k/rr);
1 = (sdf(+1))*(u(+1)*rrk(+1) + 
     q(+1)*((1 - deltau(+1)) - (phik/2)*(inv(+1)/k - delta0)^(2) + phik*(inv(+1)/k - delta0)*(inv(+1)/k)))/q;
 1/q = 1 - phik*(inv/k(-1) - delta0);
// inv = invss;
phip*(pie/piess - 1)*(pie/piess) = (1 - thetamu) + thetamu/mu + 
(sdf(+1))*phip*(pie(+1)/piess - 1)*(y(+1)/y)*(pie(+1)/piess);
// log(pie/piess) = beta*log(pie(+1)/piess)+(thetamu-1)/phip*log((1/mu)/(1/muss));
profit = (mu - 1)*y - fixedcost;
expre = (de(+1) + pe(+1))/pe;
expre2 = (de(+1) + pe(+1))^(2)/pe^(2);
varexpre = expre2 - (expre)^(2);
a = (1 - rhoa)*ass + rhoa*a(-1) + vola(-1)*ea;
vola = rhovola*vola(-1) + (1 - rhovola)*volass + volvola*evola;
z = (1 - rhoz)*zss + rhoz*z(-1) + volzss*ez;

// dc = (c^(eta)*(1 - n)^(1 - eta)-c(-1)^(eta)*(1 - n(-1))^(1 - eta))/(c(-1)^(eta)*(1 - n(-1))^(1 - eta));
// cre  = dc(+1);
// cre2 = dc(+1)^2;
// varcre = cre2 - (cre)^2;
// cre3 = dc(+1)^3;
// skewcre = cre3-3*cre2*cre+2*cre^3;
// cre4 = cre^4;
// kurtcre = cre4 - 4*cre3*cre - 3*(cre2)^2 + 12*cre2*cre^2-6*cre^4;
// mpc1 = kurtcre(+1);
// mpc2 = mpc1(+1);
// mpc3 = mpc2(+1);
// mpc4 = mpc3(+1);
// mpc5 = mpc4(+1);
// mpc6 = mpc5(+1);
// mpc7 = mpc6(+1);
// mpc8 = mpc7(+1);

// marginalutility = vf(-1)^(1/ies)*beta*expvf1sigma(-1)^((1-1/ies)/(1-sigma)-1)*vf^(-sigma)*vf^(1/ies)*eta*(1-beta)*c^(eta*(1-1/ies)-1)*(1-n)^((1-eta)*(1-1/ies));
// // cn = c^(eta)*(1 - n)^(1 - eta);
// marginalutility2 = (1+c)^4;

// expQnum2_aux = (sdf)*(((de + pe)/pe(-1))^2);
// expQnum_aux = (sdf)*((de + pe)/pe(-1));
// expQnum2 = expQnum2_aux(+1);
// expQnum = expQnum_aux(+1);
// exsdf = (sdf(+1));
// risk_neutralR2 = expQnum2/exsdf;
// risk_neutralR = expQnum/exsdf;
// risk_neutral_variance  = risk_neutralR2 - risk_neutralR^2;
// vrp = risk_neutral_variance - varexpre; 

// //expvfconbundle1 = vf(+1)^(1/ies - sigma)*a(+1)^(1 - 1/ies)*c(+1)^(eta*(1 - 1/ies) - 1)*(1 - n(+1))^((1 - eta)*(1 - 1/ies));

// //expvfconbundleaux = vf(+1)^(1/ies - sigma)*a(+1)^(1 - 1/ies)*c(+1)^(eta*(1 - 1/ies) - 2)*(1 - n(+1))^((1 - eta)*(1 - 1/ies));

// //expvfconbundle2 = vf(+1)^(2/ies - 1 - sigma)*(a(+1)^(1 - 1/ies)*c(+1)^(eta*(1 - 1/ies)-1)*(1 - n(+1))^((1 - eta)*(1 - 1/ies)))^2;

// //d2vfdc2 = eta^2*beta^2*(1/ies)*vf^(2/ies - 1)*expvf1sigma^(2*(sigma - 1/ies)/(1 - sigma))*expvfconbundle1^2 + eta^2*(sigma - 1/ies)*beta*vf^(1/ies)*expvf1sigma^((2*sigma - 1/ies - 1)/(1 - sigma))*expvfconbundle1^2 + 
// //			eta^2*(1/ies - sigma)*beta*vf^(1/ies)*expvf1sigma^((sigma - 1/ies)/(1 - sigma))*expvfconbundle2 + 
// //			eta^2*(eta*(1 - 1/ies) - 1)*beta*vf^(1/ies)*expvf1sigma^((sigma - 1/ies)/(1 - sigma))*expvfconbundleaux;

//dvfdc = eta*beta*vf^(1/ies)*expvf1sigma^(sigma - 1/ies)/(1 - sigma)*expvfconbundle1;
// //dvfdc = eta*beta*vf^(1/ies)*expvf1sigma^((sigma - 1/ies)/(1 - sigma))*expvfconbundle1;



end;

initval;
a 	 		= 	ass;
c 	 		= 	css;
de 	 		= 	dess;
deltau 		= 	delta0;
deltauprime = 	delta1;
expre 	 	= 	express;
expre2 	 	= 	expre2ss;
expvf1sigma = 	expvf1sigmass;
inv 	 	= 	invss;
k 	 		= 	kss;
n 	 		= 	nss;
mu 	 		= 	muss;
pe 	 		= 	pess;
pie 		= 	piess;
profit 		= 	profitss;
q 	 		= 	qss;
r 	 		= 	rss;
rr 	 		= 	rrss;
rrk 	 	= 	rrkss;
sdf			= 	sdfss;
u			= 	uss;
varexpre 	= 	varexpress;
vf 	 		= 	vfss;
vola 	 	= 	volass;
w 	 		= 	wss;
y	 		= 	yss;
z			= 	zss;
// sdfplus = rr*sdf;
// dc = 0;

// cre  = dc;
// cre2 = dc^2;
// varcre = cre2 - (cre)^2;
// cre3 = dc^3;
// skewcre = cre3-3*cre2*cre+2*cre^3;
// cre4 = cre^4;
// kurtcre = cre4 - 4*cre3*cre - 3*(cre2)^2 + 12*cre2*cre^2-6*cre^4;
// mpc = kurtcre;
// marginalutility = vf^(1/ies)*beta*expvf1sigma^((1-1/ies)/(1-sigma)-1)*vf^(-sigma)*vf^(1/ies)*eta*(1-beta)*c^(eta*(1-1/ies)-1)*(1-n)^((1-eta)*(1-1/ies));
// // cn = c^(eta)*(1 - n)^(1 - eta);
// mpc1 = kurtcre;

// marginalutility2 = (1+c)^4;

// mpc2 = mpc1;
// mpc3 = mpc2;
// mpc4 = mpc3;
// mpc5 = mpc4;
// mpc6 = mpc5;
// mpc7 = mpc6;
// mpc8 = mpc7;

//          expQnum2_aux = (sdf)*(((de + pe)/pe)^2);
//         expQnum_aux = (sdf)*((de + pe)/pe);
//        expQnum2 = expQnum2_aux;
//  expQnum = expQnum_aux;
// exsdf = (sdf);
// risk_neutralR2 = expQnum2/exsdf;
// risk_neutralR = expQnum/exsdf;
// risk_neutral_variance  = risk_neutralR2 - risk_neutralR^2;
// vrp = risk_neutral_variance - varexpre; 
// //expvfconbundle1 = vf^(1/ies - sigma)*a^(1 - 1/ies)*c^(eta*(1 - 1/ies) - 1)*(1 - n)^((1 - eta)*(1 - 1/ies));

// //expvfconbundleaux = vf^(1/ies - sigma)*a^(1 - 1/ies)*c^(eta*(1 - 1/ies) - 2)*(1 - n)^((1 - eta)*(1 - 1/ies));

// //expvfconbundle2 = vf^(2/ies - 1 - sigma)*(a^(1 - 1/ies)*c^(eta*(1 - 1/ies)-1)*(1 - n)^((1 - eta)*(1 - 1/ies)))^2;

// //d2vfdc2 = eta^2*beta^2*(1/ies)*vf^(2/ies - 1)*expvf1sigma^(2*(sigma - 1/ies)/(1 - sigma))*expvfconbundle1^2 + 
// //			eta^2*(sigma - 1/ies)*beta*vf^(1/ies)*expvf1sigma^((2*sigma - 1/ies - 1)/(1 - sigma))*expvfconbundle1^2 + 
// //			eta^2*(1/ies - sigma)*beta*vf^(1/ies)*expvf1sigma^((sigma - 1/ies)/(1 - sigma))*expvfconbundle2 + 
// //			eta^2*(eta*(1 - 1/ies) - 1)*beta*vf^(1/ies)*expvf1sigma^((sigma - 1/ies)/(1 - sigma))*expvfconbundleaux;
// dvfdc = eta*beta*vf^(1/ies)*expvf1sigma^(sigma - 1/ies)/(1 - sigma)*expvfconbundle1;
// //dvfdc = eta*beta*vf^(1/ies)*expvf1sigma^((sigma - 1/ies)/(1 - sigma))*expvfconbundle1;

end;

steady;

check;

shocks;
var ea 		= 1;
var evola 	= 1;
var ez	 	= 1;
end;

stoch_simul(order=4,periods=0,irf=0,noprint,nograph,nomoments,nofunctions,nocorr,pruning);











